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Abstract 

Using the concept of the Boolean derivative we study damage spreading 
for one dimensional elementary cellular automata and define their maximal 
Lyapunov exponent. A random matrix approximation describes quite well 
the behavior of "chaotic" rules and predicts a directed percolation-type phase 
transition. After the introduction of a small noise elementary cellular au- 
tomata reveal the same type of transition. 
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1 Introduction 



The behavior of the distance between two configurations submitted to the same 
dynamics (damage spreading) is considered to be a good tool to investigate the 
ergodic properties of the dynamics of discrete statistical models |y. Although the 
relation between these properties and "chaotic" behavior is still unclear, there is 
an intuitive connection between "chaos" and damage spreading on one side, and 
between a periodic attractor and damage collapsing on the other. For continuous 
dynamical systems a positive maximal Lyapunov exponent (MLE) implies chaotic 
motion. The MLE is roughly defined as the rate of the exponential divergence of 
the distance between two initially close trajectories, in the limit of long times and 
vanishing initial distances. In what follows we show how Boolean derivatives may 
be used to define the MLE of a cellular automaton. 

A Boolean one-dimensional cellular automaton (CA) is a discrete dynamical 
system defined on a lattice. The state of the system is represented by a configuration 
x = (xi, . . . , Xi, . . . , Xl) of Boolean variables, where L is the size of the lattice. We 
always use periodic boundary conditions (xj+L = The time evolution of the 
system is given by a deterministic Boolean function F 

x(t + l) = F(x(i)), (1) 

which is in turn defined by a local, uniform rule / 

Xi(t "I - 1) fi^i—rt • • • i ^ii • ■ ■ i ^i+r)i (2) 

where r is the range of the function /. There are 2 2r+1 different CA of range r. 
In what follows we restrict our study to elementary CA for which r = 1, using 
Wolfram's labeling convention ||. 

In the context of CA where time, space and dynamical variables are discrete, we 
cannot extend directly the definition of Lyapunov exponents 0, f|, || ]6|, [/J ||. Due 
to the finite interaction range r and to the finite number of states of the variables 
x — i, the distance between two initially close configurations can increase at most 
linearly for long times. 

Instead of looking at the long time behavior of the distance between two configu- 
rations, we can use some hints from the theory of continuous dynamical systems and 
study the local stability of a single trajectory with respect to a small perturbation 
(a damage or defect in the configuration). This defect can be readily recovered, or 
it can freeze being replicated without change, or finally it can propagate increasing 
the distance between the configurations. 

In section 2, we show that the distance between two configurations after the 
introduction of a defect is given by the Boolean Jacobian matrix of the evolution 
function F. While in the actual evolution of the automaton the defects can interact 
and annihilate themselves, we are interested in the stability of a single trajectory, 
and we restrict to the case of non-interacting defects which is equivalent to consider 
a product of Boolean Jacobian matrices on a trajectory. For the elementary CA it 
is a Jacobi matrix with elements equal to zero or one on the three main diagonals. 
This in turn suggests a relation with the product of random matrices of the same 
type. The MLE of the product of these random matrices shows a transition related 



to that of directed percolation ITT], 12 . 
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As reported in section 3, the results of simulations of "chaotic" CA (whose space- 
time patterns, starting from a random configuration are disordered and aperiodic) 
agree quite well with the predictions of the random matrix approximation. Adding 
a noise on the evolution of the automaton, our approach reveals the existence of a 
transition in the space of CA rules from a "frozen phase" where damage does not 
spread, to a phase where damage spreads locally with a positive MLE close to the 
one given by the product of random matrices. 



2 Boolean derivatives, defects and random matri- 
ces 

We are interested in the local stability with respect to a small perturbation, of 
the time evolution ([]]) of the configuration x. Let us denote a defect at site i the 
configuration zW with elements Zj = 5ij, j — 1, . . . , L and Sij the usual Kronecker 
symbol. The configuration y(t) = x(i) © zW differs from x(t) only at site i, where 
the XOR operation (©) is performed site by site. Depending on F and on the 
configuration x, the defect zW can originate in one time step up to three defects in 
sites i — 1. i and i + 1. We write 



x(t + 1) © y(t + 1) = A z^ 1 ) © fit A *W © F' A 



where F[, = 0, 1 and the symbol A represents the AND operation, to be performed 
between the number F[- and each element of the configuration z. The quantities 

F , d Xi (t+l) 



d Xj (t) 



are the elements of the Boolean Jacobian matrix F' of F; they are expressed by the 
Boolean derivatives of the local evolution rule / of eq. @ || pL0|| . For instance 

dxi(t +1) 

K = / \ x i-li x ii x i+l ) © / \ x i-li x ii i x i+l © !)• 

OXi+l 

As / has range 1, dxi(t + l)/dxj vanishes if \i — j\ > 1, and F' is a Jacobi matrix 
whose elements are zero or one. If the local evolution rule is expressed in terms of 
AND and XOR operations (ring sum expansion), the Boolean derivatives extract 
the linear part of /. 

We are interested in the limit of a small initial perturbation to a given trajec- 
tory. This limit corresponds in discrete dynamics to the presence of only one point 
defect. If, during the evolution, m defects appear, we consider m of replicas of 
the system and assign one of the defects to each one. We indicate with Ni(t) the 
number of replicas carrying the defect z^ at time t. If for instance we start at 
time zero with only one defect at some site i (A^(0) = 1), and the rule allows the 
spreading of the defects to the sites of the neighborhood at each time step, at t — 1 
A^_i = Ni = N i+1 = 1, at t = 2 Ni-z = N l+2 = 1, = N l+1 = 2,N l = 3, etc. 

The time evolution of the number of defects at site i is given by 

N i (t+i) = J2 F !A t ) N M 
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or, in matrix form, 

N(t+l)=F'N(t), (3) 

where the elements of F' are now interpreted as integer numbers. It is worth noting 
that Ni(t) is also the number of paths in defect space that reach the site i at time t 
starting from any defect at time t — 0. 

We define the finite-time MLE A(T) of the mapping ([3]) as 

A(T) = i2>gi7(t). (4) 
1 t 

where the local expansion rate of defects rj is defined as 

r ? (t) = |N(t+l)|/|N(t)|, (5) 

and |N| = J2i -A^- In the following A(oo) will be denoted simply by A. This definition 
is meaningful because the number |N(£)| can diverge exponentially. 

If A < the number of defects (the damage) decreases exponentially to zero, 
while if A > the damage spreads. Let us give some simple examples. Rule 
that maps all the configurations to the configuration {0} L , has A = — oo because 
the Jacobian is zero. The "chaotic" rule 150 has A = log 3, because all its Boolean 
derivatives are equal to one. A marginal case is rule 204, for which F' is the identity 
matrix and A = 0. The derivatives of the 88 "minimal" elementary CA may be 
found in ref. ||. 

In the spreading reasonable approximation to the dynamics of defects 

(0) consists in substituting the deterministic matrix F' with a random matrix of 
the same form. We therefore consider the product of random tridiagonal matrices 
M(p) having a fraction p of elements on the three principal diagonals equal to one. 
The quantity p is interpreted as the geometric mean p of the derivative on the CA 
configuration for large T, i.e., 

MT)= (JW*)J 

and 

6l i= i k =-i 

The evolution of the number of defects in the random matrix approximation 
defines a directed bond percolation problem with control parameter p, assuming 
that a site at location i at time t is "wet" if iVj(t) > 0. Observe that iVj(t) gives the 
number of directed paths that reach site i at time t inside the percolating cluster. 
Therefore we expect a second-order phase transition at p = p c with order parameter 
the density of wet sites p(t). 

We have first localized the percolation threshold at p c = 0.441(1) (where the 
number in parenthesis is the error on the last significant digit) by looking at the 
asymptotic behavior of the order parameter. Then, starting with an initial condition 
where all the sites are wet, we have verified that p(t) ~ t _/3 ^ii at p c with (3/v\\ = 
0.155(3) the usual exponents of directed percolation. 
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Figure 1: The curve shows the MLE A of a random tridiagonal matrix as a function 
of p. The diamonds show the asymptotic value of A for "chaotic" CA. Results for 
CA were obtained by letting each automaton evolve during 5000 time steps on a 512 
sites lattice and ao = 0.5. 



The results of the random matrix approximation are reported in Fig. 1, where 
the curve shows the dependence of the MLE A for the product of random matrices 
M(p) as a function of p. For p < p c A = — oo. At p c , for sufficiently large T, 

A(T) = A + aT' x 

with A = 0.237(2), \ = 0.68(4) and a > 0. This shows that at the critical point 
the number of walks on the percolation cluster grows exponentially with time, with 
an effective coordination exp(A). The exponent Xi which is not usually defined 



in percolation, might be related to the critical exponents for directed walks |]T3 
The data at the percolation threshold were obtained by letting a 10 4 x 10 4 random 
tridiagonal matrix evolve during 4000 time steps for 30 realizations. 

We obtain a mean field approximation replacing M with a constant tridiagonal 
matrix with elements equal to p. The corresponding MLE is A = log3p which is 
positive for p > 1/3. This approximation agrees well with numerical simulations for 
P > Pc with a maximum deviation of 18% at p — p c . 

In the numerical calculation of the Lyapunov exponent A one needs to renor- 
malize N(£) |14|]. This is impossible if N is defined over integers. However, since 



the Lyapunov exponents are independent of the choice of the norm ]T5[, we let 
F' (or M) act on some abstract "tangent" space in R L , using the usual Euclidean 
norm. Applying standard methods one obtains the Lyapunov exponents related to 
the exponential divergence of the norm of the product of vF^F where F ;t is the 
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Hermitian conjugate of F'. The standard definition of the MLE coincides with the 
one of eq.(|]). 

3 Elementary cellular automata 

We computed the mean number of ones n{T) in the Jacobian matrix and the finite- 
time MLE A(T) for all the 88 "minimal" elementary CA for L = 256 and L = 512 and 
5000 < T < 15000 starting from random initial configurations with a fixed fraction 
«o of live sites, «o = L^ 1 J2 x i(fy- The quantities /x(T) and A(T) are generally 
already asymptotic for T ~ 5000; moreover they show a very weak dependence on 
ao for 0.2 < cto < 0.8 (only rules 6, 25, 38, 73, 134 and 154 vary more than 10% but 
less than 20%). 

We note that 

• i) CA with constant F' independent of the configuration (e.g. rules 0, 15, 60, 
90, 150, 204) have A = log 3// with [m = 0, 1/3, 2/3 or 1. 

• ii) CA for which all configurations are mapped to a homogeneous state (e.g. 
rules 8, 32, 40, 136) have A = — oo. The control parameter /i is zero. These 
are class 1 CA in Wolfram's classification [H. 

• iii) "Chaotic" class 3 CA with nonconstant F' (e.g. rules 18, 22, 30, 41, 45, 
106, 110, 122, 126, 146) have \i > p c , A > and the damage spreads. 

The values of MLE for the "chaotic" CA of cases i) and iii) the value of A agrees 
well with the random matrix approximation, as shown in Fig. 1. This is also trivially 
true for rules of case ii. 

For the automata whose evolution leads to a nonhomogeneous periodic space 
pattern (e.g. class 2 CA), A is the logarithm of the largest eigenvalue of the prod- 
uct of the Jacobian matrices over the periodic state. The measured value of A is 
always non negative. This suggests that the asymptotic state is unstable (A > 0) 
or marginally stable (A = 0). One can think that the "freezing" of the evolution 
occurs because there are no "close" configurations which can be used as an inter- 
mediate state towards a more stable state. Therefore, we "heated" the evolution by 
exchanging the states of a small number s of pairs of randomly chosen sites at each 
time step. We observe that all the rules which had /i < 1/3 (the critical point in the 
mean field approximation) relax to /i = and A = — oo, those with 1/3 < fi < p c go 
to A = and the rules with [i > p c approach the random matrix result. In Fig. 2, 
we show the location of all the 88 minimal elementary rules with noise in the (//, A) 
diagram, starting with ao = 0.5. 

After the introduction of noise, the CA rules can be divided roughly in three 
groups, according to the value of their MLE. In the first group, with A = — oo and 
fi = we find all class 1 CA (rules 0, 8, 32, 40, 128, 136, 160 and 168) and some 
class 2 CA (rules 1, 3, 5, 7, 11, 13, 14, 19, 23, 43, 50, 72, 77, 104, 142, 178, 200 and 
232). Rules 50, 77 and 178 show very long transients of the order of 15,000 time 
steps. Rule 232, a majority rule, illustrates well a typical behavior. Configurations 
{0} L and {1} L are fixed points for this rule. A single defect in these configurations is 



6 



1.2 r- 
1.1 - 



0.9 - 

0.8 - 

0.7 - 

X 0.6 - 

0.5 - 

0.4 - 

0.3 - 

0.2 - 

0.1 - 

- 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

H 

Figure 2: The curve is the same as the one shown in Fig. 1. The diamonds show 
the values of \x and A for all the minimal CA with A > 0. 

recovered in one time step. On the other hand, an arbitraty initial configuration will 
give in a few time steps a pattern of strips. By adding a noise as described above, the 
borders of the strips perform a sort of random motion, thus allowing their merging. 
Finally, one of the two fixed points is reached, according to the initial density of the 
configuration. 

In the second group we find the "chaotic" class 3 rules. The values of \i and A 
are slightly affected by the noise. They have \i > p c and are close to the curve of the 
random matrix approximation for A. We also find CA which are not class 3 (rules 
6, 9, 25, 26, 27, 28, 29, 33, 37, 38, 54, 57, 62, 73, 94, 134, 154, and 156) but show 
local damage spreading. 

The third group contains exceptions to the random matrix approximation. This 
occurs for rules with a value of \i without noise close to 1/3 (CA 15, 51, ...204) or p c 
(CA 1,3,5,11,14,43,142,24,44,46,56,74,108 and 152) Contrary to the prediction of the 
random matrix approximation N does not vanish. Moreover, it should be remarked 
that CA with a large number of conserved additive quantities (rules 3, 4, 5, 10, 12, 
15, 34, 42, 51, 76, 138, 140, 170, 200 and 204) @ have indeed fi = 1/3 and A = 0. 

In this letter we have shown how the maxumum Lyapunov exponent can be 
defined for CA using the Boolean derivative. A positive Lyapunov exponent is 
associated to local damage spreading and on the other hand reflects the exponential 
growth of paths on the percolation clusters. For CA with < fi < p c and a positive 
Lyapunov exponent the introduction of a small noise produces the collapse to A = 
or A = — oo. A random matrix model is directly suggested by the CA dynamics 
and displays a directed percolation phase transition. A somewhat similar phase 
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transition is observed in CA rule space in the presence of a small amount of noise. 
The extension of our definition of Lyapunov exponent to other discrete systems, and 
possibly to probabilistic dynamics will be the subject of future investigations. 
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